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Maximally dissipative boundary conditions are applied to 
the initial-boundary value problem for Einstein's equations in 
harmonic coordinates to show that it is well-posed for homo- 
geneous boundary data and for boundary data that is small in 
a linearized sense. The method is implemented as a nonlinear 
evolution code which satisfies convergence tests in the nonlin- 
ear regime and is stable in the weak field regime. A linearized 
version has been stably matched to a characteristic code to 
compute the gravitational waveform radiated to infinity. 

The waveform emitted in the inspiral and merger of a 
relativistic binary is theoretical input crucial to the suc- 
cess of the fledgling gravitational wave observatories. A 
computational approach is necessary to treat the highly 
nonlinear regime of a black hole or neutron star colli- 
sion. Developing this computational ability has been the 
objective of the Binary Black Hole (BBH) Grand Chal- 
lenge [1] and other world wide efforts. The Grand Chal- 
lenge based a Cauchy evolution code on the Arnowitt- 
Deser-Misner (ADM) [2] formulation of Einstein's equa- 
tions. Exponentially growing instabilities encountered 
with that code have been traced to improper bound- 
ary conditions [3]. (Even in the absence of boundaries, 
an ADM system linearized off a Minkowski metric has 
a power law instability [4]; this triggers an exponen- 
tial instability when the background Minkowski metric 
is treated in non-Cartesian, e.g. spherical, coordinates.) 
Other groups have encountered difficulties in treating 
boundaries (see [5] for a recent discussion) and the work- 
ing practice is to forestall problems by placing the outer 
boundary far from the region of physical interest (see 
e.g. [6]) or compactify the spacetime (see e.g. [7]). 

This deficiency extends beyond numerical relativity to 
a lack of analytic understanding of the initial-boundary 
value problem (IBVP) for general relativity. The local- 
version of the IBVP is schematically represented in Fig. 1. 
Given Cauchy data on a spacelike hypersurface S and 
boundary data on a timelike hypersurface B, the problem 
is to determine a solution in the appropriate domain of 
dependence. Whereas there is considerable mathematical 
understanding of the gravitational initial value problem 
(for reviews see [8]), until recently the IBVP has received 
little attention (see e.g. [9,10]). Indeed, only relatively re- 
cently has the method of maximally dissipative boundary 
conditions [11,12] been extended to the nonlinear IBVP 
with boundaries containing characteristics [13,14], such 
as occurs in symmetric hyperbolic formulations of gen- 



eral relativity. Friedrich and Nagy [10] have applied these 
methods to give the first demonstration of a well-posed 
IBVP for Einstein's equations. The Friedrich-Nagy work 
is of seminal importance for introducing the maximally 
dissipative technique into general relativity. Their formu- 
lation, which uses an orthonormal tetrad, the connection 
and the curvature tensor as evolution variables, is quite 
different from the metric formulations implemented in 
current codes designed to tackle the BBH problem. Al- 
though it is not apparent how to apply the details of 
the Friedrich-Nagy work to other formalisms, the general 
principles can be carried over provided Einstein's equa- 
tions are formulated in the symmetric hyperbolic form 



J2A a {u)d a u = S(u) 



(0.1) 



with coordinates x a = (t,x,y,z) = {t,x l ) and evolution 
variables u = (iti, ujy), where A a are NxN symmetric 
matrices and A* is positive-definite. 

The simplest symmetric hyperbolic version of Ein- 
stein's equations employs harmonic coordinates satisfy- 
ing H a := y/—g\2x a = df){y/—gg a P) = 0, in which the 
well-posedness of the initial value problem was first es- 
tablished [15,16]. Well-posedness expresses the existence 
of a unique solution with continuous dependence on the 
data. In the nonlinear case, existence is only guaranteed 
for a short time, reflecting the possibility of singularity 
formation. Here we show how this approach (i) can be 
applied to the IBVP problem in harmonic coordinates, 
(ii) can be implemented as a robustly stable, convergent 
3-dimensional nonlinear Cauchy evolution code and (iii) 
can be accurately matched, in the linearized approxi- 
mation, to an exterior characteristic evolution code to 
provide the proper physical boundary condition for com- 
puting the waveform radiated to infinity by an isolated 
source. Reference [17] discusses the suitability of har- 
monic coordinates for numerical work and for simulating 
the approach to a curvature singularity. 

We base the evolution on the metric density 7 a/3 = 

v"/3 



with g = det(7°^) = det(g a p). As in 



Eq. (B.87) of Fock [18], we split the Einstein ten- 
sor into G af} = £ af3 + \g a0 B - B aP , where B a(3 = 
— (—g)^ 1 / 2 y^ a H^ vanishes when H a = 0, where 
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and where S a ^ contains no second derivatives of the 
metric. When the harmonic conditions H a = are 
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satisfied, the reduced Einstein equations £ a ^ = have 
principal part which is governed by the wave opera- 
tor ^i^d^dy. In terms of the evolution variables u — 
^ a f3 ^ T ap ^ x a P ,y<*P ,Z a P), where T af3 = da aP , X a/3 = 
dxT'\ y afi = dyT li and Z a(i = <9 z7 a/3 , we put these 
wave equations in symmetric hyperbolic form (0.1) by a 
standard construction [19]. 

We adapt the harmonic coordinates to the boundary 
so that the evolution region lies in z < 0, with the bound- 
ary fixed at z — in the numerical grid. We write 
x a = (x a ,z) = (t,x,y,z) to denote coordinates adapted 
to the boundary, so that x a = (x a ,z) = (t,x l ) depend- 
ing whether the Latin index is near the beginning or end 
of the alphabet. We further adapt the coordinates so 
that "f za \i3 = (and hence T za \ B = 0). For any timclikc 
B, these harmonic gauge conditions can be satisfied and 
they are assumed throughout the following discussion. 

We base the well-posedness of the homogeneous IB VP 
on Theorems 2.1 and 2.2 of Secchi [14], which require 
that u satisfy a boundary condition of the form Mu = 0, 
where M is a matrix independent of u and of maximal 
rank such that the normal flux T z = (u, A z u) associated 
with an energy norm be non-negative. Secchi's theorems 
include the present case where the boundary is "charac- 
teristic with constant rank" , i.e. A z has a fixed number 
of eigenvalues. The above gauge conditions consider- 
ably simplify this maximally dissipative condition for the 
reduced system £ a/3 = 0. The boundary matrix takes the 
form A z = -f zz C, where C is a constant matrix, and the 
flux inequality reduces to 

T z = -^T at) Z afi > 0. (0.3) 

Here T z is identical to the standard energy flux for 
the sum of 10 independent scalar fields. This require- 
ment can be satisfied in many ways, e.g. by combina- 
tions of the homogeneous Dirichlet boundary condition 
dtl — T a P — 0, the homogeneous Neumann condition 
(9 z <-y a/3 = Z a/3 = and the homogeneous Sommcrfcld con- 
dition (d t + d z y f al3 = T af) + Z a(i = on the various field 
components. All these boundary conditions have the re- 
quired form Mu = 0. The maximality of the rank of M 
ensures that boundary conditions only be applied to vari- 
ables propagating along characteristics entering the evo- 
lution region from the exterior [20]. For instance, assign- 
ment of a boundary condition to the variable T Q/3 — Z al3 , 
which propagates from the interior toward the boundary, 
would violate (0.3). 

Whereas the well-posedness of the IBVP for the re- 
duced system can be accomplished by a variety of bound- 
ary conditions, it can only be established for the full sys- 
tem in a limited sense. The Bianchi identities and re- 
duced equations imply V p (B All/ — ^g^B) = 0, which has 
the explicit form 

^ v d^d v H a + C^d^H" + D"H V = 0, (0.4) 

where and I?" depend algebraically on u and H a . 
Thus H a obeys a symmetric hyperbolic equation of the 



form (0.1). Harmonic Cauchy data satisfying the Hamil- 
tonian and momentum constraints and H a = also sat- 
isfy dtH a — on S by virtue of the reduced equations, so 
that uniqueness guarantees H a = in the domain of de- 
pendence T>i (see Fig. 1) and hence the well-posedness of 
the Cauchy problem for the full system. To extend well- 
posedness to the homogeneous IBVP, i.e. to include re- 
gion T>2, we impose boundary conditions for the reduced 
system that imply the homogeneous boundary conditions 
H z = and d z H a — for the harmonic constraints. 
Combined with the gauge condition j za = 0, the con- 
dition H z := 9h7 zb + d z Y' z = requires the Neumann 
boundary condition Z zz = 0. We also impose the homo- 
geneous Neumann boundary conditions Z ab = so that 
d z H a := d b Z ab + d 2 z -i az = requires d^ az = at the 
boundary. Remarkably, subject to the above conditions, 
the reduced equation £ az = implies d z j az — at the 
boundary! Underlying this result is that S az — at the 
boundary due to the local reflection symmetry implied 
by the above conditions. This establishes the maximally 
dissipative boundary conditions H z = d z H a = for the 
constraint propagation equations (0.4) which ensure that 
the full Einstein system is satisfied. 

In practice, homogeneous boundary conditions do not 
correspond to a given physical problem, e.g. homoge- 
neous Neumann data at the end of a string lead to a 
free endpoint whereas the the endpoint might be under- 
going a forced oscillation requiring inhomogeneous data. 
This flexibility is supplied within the maximally dissi- 
pative formalism by the ability to extend the homoge- 
neous boundary condition Mu = to the inhomogeneous 
form M(x a )(u — q(x a )) = [10]. This preserves the well- 
posedness of the IBVP for the reduced system with in- 
homogeneous Neumann data Z zz — q zz and Z ab = q ab . 
For the full system, the gauge condition j za = and 
the boundary constraint H z = forces q zz — 0. Next, 
d z H a = implies 

D a (Z ab /^g~^Tg) = 0, (0.5) 

where D a is the connection intrinsic to the boundary. 
The appearance of the metric and D a in Eq. (0.5) intro- 
duces M-dependence in the boundary data so that Sec- 
chi's theorems do not apply. However, the theory does 
apply to boundary data Sq ab linearized off a nonlinear so- 
lution with homogeneous data, either exact or generated 
numerically. Then Eq. (0.5) has the form 

d b Sq ab + F« c (x d )5q bc -0, (0.6) 

where F b a c (x d ) is explicitly known via the metric and 
connection of the homogeneous solution. The principal 
part of Eq. (0.6) is identical to the analogous equation 
in the linearized version of the harmonic IBVP treated 
in Ref. [21]. In terms of the coordinates x a — (t,x A ) — 
(t, x,y) on the boundary, a simple transformation of vari- 
ables (see [21]) recasts Eq. (0.6) as a symmetric hy- 
perbolic system of equations for (ft — \5AB^q AB and 
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y A = Sq tA . (Here Sab is the Kronecker delta.) This 
system uniquely determines <f> and y A in terms of their 
initial values and the 3 pieces of free boundary data 
y AB = gq AB + 5 AB( 5q tt _ S CD 5q CD ). (This in accord 

with Ref. [10], although there is no direct correspondence 
with the three free pieces of boundary data in [10]). Since 
only coordinate conditions have been imposed here, the 
only restriction on physical generality is the linearity of 
the boundary data. 

The IBVP also requires consistency between the 
Cauchy data and boundary data at S U B, which deter- 
mines the degree of differentiability of the solution [14] . 
As in the string example, consistent homogeneous Neu- 
mann boundary data and Cauchy data imply a virtual 
reflection symmetry across the boundary, which is bro- 
ken in the inhomogeneous case. Although the IBVP 
is well-posed for the reduced system and for the con- 
strained system with boundary data linearized about the 
homogeneous case, no available theorems guarantee well- 
poscdness for the constrained inhomogeneous case. In 
this respect, the analytic underpinnings are not as gen- 
eral as the Friedrich-Nagy formulation. Numerical simu- 
lations are necessary to shed further light on this ques- 
tion. The key feature of our formulation is that if a so- 
lution exists, as provided by a convergent numerical sim- 
ulation, then it necessarily satisfies the constraints, since 
the constraint propagation equation (0.4) is then satis- 
fied with maximally dissipative, homogeneous boundary 
data. In the strong field convergence tests described be- 
low, exact solutions provide the Cauchy and boundary 
data. 

In constructing a code to demonstrate these results, 
we take considerable liberty with the symmetric hyper- 
bolic formalism. In particular, we use the second dif- 
ferential order form of the equations based upon the 10 
variables j a/3 rather than the 50 first order variables it; 
we use a cubic boundary aligned with Cartesian coordi- 
nates, although the mathematical theorems only apply to 
smooth boundaries; and we replace the gauge condition 
Y' a = by 7 za = q a (x b )Y' z , where q a — d z x a \ B is the 
free Neumann boundary data in the transformation to a 
general harmonic coordinate system satisfying Ox a = 0. 
The harmonic boundary constraint H z = now implies 
q zz = —d a (q a Y~ z )\i3 and the constraint d z H a = again 
determines <j) an d y A in terms of the free boundary data 
(q a ,y AB ), now through a symmetric hyperbolic system 
obtained from adding source terms arising from q a ^ 
to the right hand side of Eq. (0.5). We use the finite 
difference techniques described in [21], where robust sta- 
bility and convergence of a linearized harmonic code was 
demonstrated. In the linearized theory, the decoupling 
of the metric components gives more flexibility in formu- 
lating a well-posed IBVP. The linearized harmonic code 
could be consistently implemented with Dirichlet bound- 
ary conditions, in which case it ran stably for 2000 cross- 
ing times even with a piecewise-cubic spherical bound- 
ary cut out of the Cartesian grid. However, we have not 
found a well-posed version of the nonlinear theory that 



avoids Neumann boundary conditions and the associated 
numerical complications which we describe below. 

We tested robust stability [23] of the nonlinear code 
by initializing the evolution with random, constraint vi- 
olating initial data 7 a/3 = -n al3 + e Qf9 and by assign- 
ing random boundary data q{x a ) = e at each point 
of the cubic grid boundary, with the e's random num- 
bers in the range (— 10~ 10 , 10~ 10 ). (Although differing 
from the standard numerical definition of stability re- 
lated to convergence, robust stability is computation- 
ally practical for revealing short wavelength instabili- 
ties.) Under these conditions, the noise in the nonlinear 
code grows linearly at the same rate for 2000 crossing 
times for both the 48 3 and 72 3 grids. We tested conver- 
gence in the nonlinear regime using a gauge-wave gener- 
ated by the harmonic coordinate transformation (a;, y) = 
x B — > x B + a B sin [2n (\/3t + x + y + z)] acting on the 
Minkowski metric, with a x = 0.06 A, a y = 0.04 A. The 
resulting gauge-wave has amplitude \ \g a(3 — rj al3 \\ 00 w A. 
We use periodicity in the (x, y) plane to evolve with 
smooth toroidal boundaries at z — ±1/2. Second-order 
convergence in the non-linear regime was confirmed with 
the amplitude A = 10 _1 . Figures 2 demonstrates the 
convergence of the solution and Fig. 3 shows the ab- 
sence of anomalous boundary error. Error arising from 
the application of Neumann boundary conditions even- 
tually triggers a nonlinear instability, which occurs after 
30 crossing times with the 120 3 grid. Runs with the am- 
plitude A = 10~ 3 were carried out on the 80 3 grid for 
300 crossing times without encountering the above non- 
linear instability (see Figure 2). In the case of a cubic 
boundary, the nonlinear code cleanly propagates a phys- 
ical pulse with amplitude 10~ 7 that corresponds to an 
exact linearized solution; but, for a gauge-wave of ampli- 
tude A = 10~ 3 , substantial error arises at the edges and 
corners due to our present method of applying Neumann 
boundary conditions and leads to an instability after 60 
crossing times. 

The physically proper boundary data for a given 
problem is a separate and difficult problem for non- 
linear systems. One approach is to supply q(x a ) by 
Cauchy-characteristic matching (CCM) in which an in- 
terior Cauchy evolution with cubic boundary is matched 
to an exterior characteristic evolution on a sequence of 
outgoing null cones extending to infinity (for a review see 
[22]). In simulations of a nonlinear scalar wave with peri- 
odic source, CCM was demonstrated to compute the ra- 
diated waveform more efficiently and accurately than ex- 
isting artificial boundary conditions on a large but finite 
boundary. [24] Previous attempts at CCM in the grav- 
itational case were plagued by boundary induced insta- 
bilities growing on a scale of 10 to 20 grid crossing times. 
Although stable behavior of the Cauchy boundary is only 
a necessary but not a sufficient condition for CCM, tests 
carried out with a linearized harmonic Cauchy code with 
a well-posed IBVP matched to a linearized characteristic 
code show no instabilities. 
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In the tests of CCM, the linearized Cauchy code was 
supplied outer boundary data q in Sommcrfcld form by 
the exterior characteristic evolution and boundary data 
for the characteristic code was supplied on an interior 
spherical boundary by the Cauchy evolution. Robust sta- 
bility for 2000 crossing times on a Cauchy grid of 45 3 was 
confirmed. For a linearized wave pulse, Figure 4 shows 
a sequence of profiles of the metric component j xy prop- 
agating cleanly through the spherical boundary as the 
wave pass to the characteristic grid, where it is propa- 
gated to infinity. Further details and tests of CCM and 
the question of its extension to the nonlinear theory will 
be reported elsewhere. 

At present, the major limitation in the nonlinear code 
stems from the difficulty in handling large values of 7* z 
at the boundary. This is evidenced by numerical exper- 
iments with the manifestly well-posed IBVP consisting 
of a scalar wave propagating between smooth toroidal 
boundaries according to the flat-space wave equation 

( - d 2 t - 2vd t d z + d 2 x + d 2 y + (i - v 2 )d 2 ^ $ = o, 

which arises from the transformation z —* z + vt on stan- 
dard incrtial coordinates. The value of 7* z represents the 
velocity of the boundary relative to observers at rest with 
respect to the Cauchy slicing. For the flat space wave 
equation in second order form, there have apparently 
been no studies of numerical algorithms which apply Neu- 
mann boundary conditions to such moving boundaries. 
In fact, only very recently has there been a thorough 
treatment of Neumann boundary conditions for the the 
flat space wave equation with a stationary (but curvilin- 
ear) boundary [25]. This treatment uses Neumann data 
to update the field at a boundary point at the current 
time step by a one-sided finite difference approximation 
for the normal derivative. Such stencils for approximat- 
ing normal derivatives apply only when the normal direc- 
tion is tangent to the Cauchy slicing, i.e. when g tz \e = 0. 
The general case in which g tz \s ^ requires a more com- 
plicated stencil involving interior points to the future or 
past of the current time step. We have developed a new 
approach which successfully handles this general case for 
the above scalar wave test problem but requires further 
refinement to handle boundaries with edges and corners 
before it can be implemented in the gravitational code. 

We thank H. Friedrich and B. Schmidt for educating 
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with help from the Cactus development team of the AEI. 
The work was supported by NSF grant PHY 9988663. 



Introduction to Current Research, ed. L. Witten (Wiley, 
NY, 1962). 

[3] B. Szilagyi, R. Gomez, N. T. Bishop and J. Winicour, 

Phys. Rev. D, 62, 104006 (2000). 
[4] G. Calabrese, J. Pullin, O. Sarbach and M. Tiglio, Phys. 

Rev. D, 66, 041501 (2002). 
[5] G. Calabrese, L. Lehner and M. Tiglio, Phys. Rev. D, 65, 

104031 (2002). 

[6] R. L. Marsa and M. W. Choptuik, Phys. Rev. D, 54, 4929 
(1996). 

[7] S. Husa, Problems and Successes in the Numerical Ap- 
proach to the Conforrnal Field Equations, in The Confor- 
mal Structure of Spacetimes: Geometry, Analysis, Nu- 
merics, ed. J. Frauendiener and H. Friedrich, Lecture 
Notes in Physics 604 (Springer, New York, 2002). 

[8] O. Reula, Living. Rev. ReL, 1, 3 (1998); H. Friedrich and 
A. Rendall, in Einstein's Equations and Their Physical 
Implications, edited by B. Schmidt (Springer, New York, 
2000); A, Rendall, Living. Rev. ReL, 4, 1 (2001). 

[9] J. M.Stewart, Class. Quantum Grav., 15, 2865 (1998). 
[10] H. Friedrich and G. Nagy, Commun. Math. Phys., 201, 
619 (1999). 

[11] K. O. Friedrichs, Comm. Pure Appl. Math., 11, 333 
(1958). 

[12] P. D. Lax and R. S. Phillips, Comm. Pure Appl. Math., 

13, 427 (1960). 
[13] J. Rauch, Trans. Am. Math. Soc, 291, 167 (1985). 
[14] P. Secchi, Arch. Rational Mech. Anal., 134, 155 (1996). 
[15] Y. Foures-Bruhat, Acta. Math., 88, 141 (1955). 
[16] A. E. Fisher and J. E. Marsden, Comm. Math. Phys., 28, 

1 (1972). 

[17] D. Garfmkle, Phys.Rev., D65, 044029 (2002). 

[18] V. Fock, The Theory of Space, Time and Gravitation, p. 
392 (MacMillan, New York, 1964). 

[19] R. Courant and D. Hilbert, Methods of Mathematical 
Physics, Vol. II, p. 594 (Interscience, NY, 1962). 

[20] B. Gustafsson, H.-O. Kreiss and J. Oliger, Time De- 
pendent Problems and Difference Methods (Wiley, NY, 
1995). 

[21] B. Szilagyi, B. Schmidt and J. Winicour, Phys. Rev., 

D65, 064015 (2002). 
[22] J. Winicour, Living. Rev. ReL, 4, 3 (2001). 
[23] The robust stability test has been extended to include 

dependence on grid size in "Toward standard testbeds for 

numerical relativity", M. Alcubierre et al, gr-qc/0305023. 
[24] N. T. Bishop et al, J. Comp. Phys., 136, 140 (1997). 
[25] H.-O. Kreiss, N. A. Petersson and J. Ystrom, SI AM 

Journal on Numerical Analysis, 40, No. 5, pp 1940-1967 

(2002). 



[1] S. Brandt et al, Phys. Rev. Lett, 85, 5496 (2000). 

[2] R. Arnowitt, S. Deser and C. Misner, in Gravitation: An 



4 




FIG. 1. Schematic representation of the domain of depen- 
dence X>i of the initial value problem and the domain of de- 
pendence T>i U V 2 of the IBVP. 



FIG. 3. A y = slice of the metric component 7" , evolved 
for 30 crossing times, amplitude A = 10 _1 , with a toroidal 
boundary in the (x, y) plane. 
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FIG. 2. The Loo norm of the finite-difference error 
ll z = "/ana ~ ittm, rescaled by a factor of 1/A 2 , for a 
gauge- wave. The upper two (mostly overlapping) curves 
demonstrate convergence to the analytic solution for a wave 
with amplitude A = 10 _1 with gridsizes 80 3 and 120 3 . We 
also plot l-ffloo, the L x norm of ^/{H 1 ) 2 + 5^WW, to 
demonstrate that convergence of the harmonic constraints is 
enforced by the boundary conditions. The lower curve rep- 
resents evolution of the same gauge-wave with A = 10~ 3 for 
300 crossing times with gridsize 80 3 . 
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FIG. 4. Sequence of z = slices of the metric component 
y xy , evolved for one crossing time, with the linear matched 
Cauchy-characteristic code. In the last snapshot, the wave 
has propagated cleanly onto the characteristic grid. 
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